#ifndef __DEME_H__
#define __DEME_H__

//#include <classlib\arrays.h>
#include "arrays.h"
#include "public.h"
#include "migrmat.h"
#include "genealgy.h"
#include "genstat.h"

#define MAXLONG	0x7fffffff

const long MAX_GENERATIONS=10000000L;
const long MIN_SIZE=1;
const long MAX_SIZE=2000000000L;

//------------------------------------------------------------------------------
//Christian 06/28/04: The TSamp class. If samples are to be taken from
//the same deme at multiple time points, sample size cannot be treated as a
//simple <long>
class TSamp {
	friend ostream& operator<<(ostream& os, const TSamp& S);
	friend istream& operator>>(istream& is, TSamp& S);

    private :
      long _size;  //number of samples
      long _age;   //number of generations old these samples are
      int _deme;   //the deme this sample comes from
      int _stat_grp; //the stat group of the sample

    public :
        TSamp () {;};
        void SetVals(const long& sz, const long& ag, int d, int sg) {
            _size=sz; _age=ag; _deme=d; _stat_grp=sg;
        };
        ~TSamp() {;};
        const long& n(void) const { return (_size); };
        const long& gen(void) const { return (_age); };
        const int& deme(void) const { return (_deme); };
        const int& stat_grp(void) const { return (_stat_grp); };
        bool operator< (TSamp t) {
            if (_age<t._age) return 1;
            return 0;
        }
};

extern ostream& operator<<(ostream& os, const TSamp& S); 
extern istream& operator>>(istream& is, TSamp& S);

typedef MY_TSArrayAsVector<TSamp>				TSampArray;  //Christian 6/29/04

//------------------------------------------------------------------------------
class TDeme {

	friend ostream& operator<<(ostream& os , const TDeme& D);

	private :
      static TMigrationMatrix* Migrations; //A pointer to a migration matrix
      float _growth_rate;
      double  _deme_size;
      long  _min_deme_size, _max_deme_size;
      TSampArray _samples;  //Christian 6/29/04
      float _proba_coal;
      long _lidum;    			//dummy variable for drawing random numbers
      int _id;

      bool _is_migration;
      long _polym_sites;
   public :
      long  _lineages;
      long _oldest_samp;
      TListNode OriginalNodes;  //A list of the nodes originally present at generation zero in the deme
      TListNode NodeList;    	//A list of the nodes currently present in the deme

   	TDeme() : NodeList(), OriginalNodes() {
         initialize();
      };
      TDeme(const TDeme& D)  : NodeList(), OriginalNodes() {
      	if (this!=&D) *this=D;
      };
   	TDeme(const long& ds, const TSampArray samps, const float& gr) : NodeList(), OriginalNodes()  {
         _growth_rate=gr;
         _samples=samps;
         _deme_size=(double)ds;
         _lineages=_samples.count();
         _lidum=1L;
         _is_migration=true; 
         Migrations=NULL;
       	_polym_sites=0L;
       	_oldest_samp=0L;
      };
      ~TDeme() {;};
      void set_growth_rate(const float& gr) {_growth_rate=gr;}
      //Christian 6/29/04: Modified for ancient sample compatibility
      /*void set_sample_size(const TSamp& s) {
      	 _sample=s;
         _lineages=s.n();
      }*/
      void AddSample(const TSamp& s) {
      	 _samples.Add(s);
         _lineages+=s.n();
      }
      void set_size_bounds(const long& min, const long& max) {
			_min_deme_size=min;
			_max_deme_size=max;
      }
      void set_deme_size(const long& size) {
      	_deme_size=(double)size;
   		if (size) _proba_coal=1.0/_deme_size;
         else _proba_coal=0.0;
      }
      void set_id(const int id) {_id=id;}
      void initialize() {
      	 _growth_rate=0;  //Stationary population
         //_samples=NULL; 
         // Christian 6/29/04: the array is should already by properly
         // initialized by MY_TSArrayAsVector
         _deme_size=1.0;
         _min_deme_size=MIN_SIZE;
         _max_deme_size=MAX_SIZE;
         _lineages=0;
         _lidum=1L;        
         _is_migration=true;
         Migrations=NULL; 
       	_polym_sites=0L;
       	_oldest_samp=0L;
      }
      void set_migration_matrix(TMigrationMatrix* const MM) {
      	Migrations=MM;
      }
      TDeme& operator=(const TDeme& D);
      int operator==(const TDeme& D) {
      	if (_id==D._id) return 1;
         return 0;
      };
      int operator<(const TDeme& D) { return _id<D._id; }
      int copy_node_pointers(const TTree& T, const int curr_lin, const long& max_lin);
      //Implementation of a single migration event towards another deme (sink)
      int migrate(TDeme& sink, int x) {
   		//Add selected node to sink and then increments sink's number of lineages
   		sink.NodeList[sink._lineages++]=NodeList[x];
   		//Decrements last node's position, and replace selected node by last node
   		//CHRISTIAN: BUG!!! What if the migrating node *IS* the last node?
   		if(x!=(--_lineages)) NodeList[x]=NodeList[_lineages];
         //printf("\n-Moved #%d from %d to %d",x,_id,sink._id);
         return 1;
      };

      //Check and implements migrations
      int send_migrants(TDeme& sink, const int& i,const long& time);

      //Implements the coalescence process within population
      int coalesce_lineages(const long& time, TTree& GeneTree, int pos, bool hudson);
      float growth() {return (_growth_rate!=0.0);}

      //Implements an exponential growth of population for positive _growth_rate
      //and and exponential decline for negative _growth_rate values
      float exponential_resize() {
         _deme_size*=exp(_growth_rate);
         if ((long)_deme_size<_min_deme_size) {
         	_deme_size=_min_deme_size;
			   _growth_rate=0.0;
         }
         if ((long)_deme_size>_max_deme_size) {
         	_deme_size=_max_deme_size;
			   _growth_rate=0.0;
         }
   		if (_deme_size) _proba_coal=1.0/_deme_size;
         else _proba_coal=0.0;
         return _deme_size;
      }

      //Implements a resizing of the deme size according to a given factor
      void linear_resize(const float& factor) {
      	_deme_size*=factor;
   		if (_deme_size) _proba_coal=1.0/_deme_size;
         else _proba_coal=0.0;
      }
      void reset();
      long sample_size() {return _samples.count();}
      int  check_migrations();
      bool is_migration() {return _is_migration;}
      long deme_size() {return (long)_deme_size;}
      void print_haplotypes_for_Arlequin(ostream & os, const int& n, const Mut_Type data_type);
      void print_haplotypes_for_PAUP(ostream & os, const Mut_Type data_type);
};

extern ostream& operator<<(ostream& os , const TDeme& D);


//------------------------------------------------------------------------------
// A class specifying the conditions of a demographic historical change
// sometimes in the past, like the fusion of two demes, or a change in
// deme size

class THistoricalEvent {
	friend ostream& operator<<(ostream& os, const THistoricalEvent& HE);
	friend istream& operator>>(istream& is, THistoricalEvent& HE);
	public:
	   long time;     			//time in the past at which the event occurs (in generations)
      int source;    			//deme which sends of migrants
      int sink;      			//deme which receives migrants
      float migrants; 			//Proportion of migrants sent to sink from source
      float new_deme_size; 	    //New relative deme size for SOURCE deme
      float new_growth_rate; 	//New growth rate for SOURCE deme
      int MigMat;				//New migration matrix number
      THistoricalEvent() {};
      ~THistoricalEvent() {};
      int operator==(const THistoricalEvent& HE) {
      	if (source==HE.source &&
     		sink==HE.sink &&
            time==HE.time &&
            migrants==HE.migrants &&
            new_deme_size==HE.new_deme_size &&
            new_growth_rate==HE.new_growth_rate &&
            MigMat==HE.MigMat) return 1;
         return 0;
      }
      int operator<(const THistoricalEvent& HE) const {
      	return time<HE.time;
      }
};

extern ostream& operator<<(ostream& os, const THistoricalEvent& HE); 
extern istream& operator>>(istream& is, THistoricalEvent& HE);

//------------------------------------------------------------------------------
typedef MY_TArrayAsVector<TDeme> 			    TDemeArray;
typedef MY_TSArrayAsVector<THistoricalEvent>    TEventArray;
typedef MY_TArrayAsVector<long>                 TDemeSizeArray;
typedef MY_TArrayAsVector<float>                TGrowthArray;
typedef MY_TArrayAsVector<TMigrationMatrix>     TMigrMatArray; //Loro_16_9_99

class TDemeCollection {
	friend ostream& operator<<(ostream& os,const TDemeCollection& DC);
  	private :
   	//TDemeArray  Demes;
      TDeme          *Demes;
      TEventArray 	Events;
      TMigrMatArray 	*MigrMatArray;
      TMigrationMatrix 	*CoalTimes; //A matrix of mean pairwise coalescence times
      TMigrationMatrix 	*PairDiff; //A matrix of mean number of pairwise differences
      TMigrationMatrix 	*MinCoalTimes; //A matrix of minimum pairwise coalescence times
      TTree			   GeneTree;
      bool _is_migration;
      int _num_demes; 
      long _polym_sites;
      long _tot_mut;
      int _cur_mig_mat;
      int _num_stat_grps; //Christian
   public :
   	TDemeCollection() : Events(1,10), GeneTree() {
      	CoalTimes=NULL;
         MinCoalTimes=NULL;
         PairDiff=NULL;
         Demes=NULL;
         _num_demes=0; 
       	 _polym_sites=0L;
       	 _tot_mut=0L;
         _cur_mig_mat=0;
         _num_stat_grps=0;
      };
      ~TDemeCollection() {
      	if (CoalTimes) delete CoalTimes;
      	if (PairDiff) delete PairDiff;
      	if (MinCoalTimes) delete MinCoalTimes;
         if (Demes) delete[] Demes;
      }
      TDemeCollection& operator=(const TDemeCollection& DC);
      long count_lineages();
      int reajust_deme_size(const int i, const long& ds);
      int stat_grps() {return _num_stat_grps; }; //Christian
      int create_demes(const TMigrMatArray* const MA);
      int create_lineages(const TSampArray& SS);
      int initialize_growth_rates(const TGrowthArray& GR);
      int initialize_deme_size(const TDemeSizeArray& DS);
      int initialize_events(const TEventArray& HE);
      int initialize_samp_size(const TSampArray& SS); //Christian 6/29/04
      int build_tree();
      TMigrationMatrix& migr_matrix() {return (*MigrMatArray)[_cur_mig_mat];}
      MatElemType& min_coal_time(const int i, const int j) const {
      	return (*MinCoalTimes)(i,j);
      }
      MatElemType& mean_pair_diff(const int i, const int j) const {
      	return (*PairDiff)(i,j);
      }
      void print_gene_tree(ostream& os, const tree_output_type tree_type, const float & mu) {
      	GeneTree.print_nodes(os, tree_type, mu);
      }
      void print_sequences(ostream& os, const Mut_Type& mut_type) {
      	GeneTree.print_sequences(os, mut_type);
      };
      int implement_event(const int& cur_event);
      void reset(const TDemeSizeArray& SA, const TGrowthArray& GA);
      double mean_coal_time() {return GeneTree.mean_coalescence_time();}
      double sd_coal_time() {return GeneTree.sd_coalescence_time();}

      const TNode& tree_node(const int& pos) {return GeneTree[pos];}
      TDeme& operator[](const int& d) {return Demes[d];}
      int print_coal_time_matrix(ostream& os ) {
      	if (!CoalTimes) return 0;
         os << "\n" << *CoalTimes;
         return 1;
      }
      int print_mean_pair_diff_matrix(ostream& os ) {
      	if (!PairDiff) return 0;
         os << "\n" << *PairDiff;
         return 1;
      }
      const TMigrationMatrix& coal_time_mat() const {return *CoalTimes;}
      const TMigrationMatrix& mean_pair_diff_mat() const {return *PairDiff;}  
      bool is_migration() {return _is_migration;}
      int check_migrations();
      int check_growth();
      int num_demes() const {return _num_demes;} 
      // Modified declaration to add 12 parameters by Rumen Kostadinov
      const long& sprinkle_mutations(const float& mut_rate, const float& new_mut_rate, const long& time_of_mut_rate_change, 
				     const int& len,
				     const Mut_Type& mut_type, const float& gamma_par,
				     const float& mut_ratio, const float& prop_sites,
				     const float& trans_rate, 
				     const float& a_c, const float& a_g, const float& a_t,
				     const float& c_a, const float& c_g, const float& c_t,
				     const float& g_a, const float& g_c, const float& g_t,
				     const float& t_a, const float& t_c, const float& t_g,
				     const float& alpha, const float& beta,
				     const int& range_constr);
      void write_samples_to_Arlequin_file(ostream& os, const Mut_Type data_type);
      void write_group_section_to_Arlequin_file(ostream& os);
      void write_samples_to_PAUP_file(ostream& os, const Mut_Type data_type);
      void write_Statfile(ostream& os, const Mut_Type data_type);
      int migr_mat() {return _cur_mig_mat;}
};

extern ostream& operator<<(ostream& os,const TDemeCollection& DC);

//------------------------------------------------------------------------------

#endif
